{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mouse-to-Human alignment/mutation modeling\n",
    "\n",
    "This notebook checks whether mutations in MSK IMPACT dataset can be modelled in mouse based on amino acid conservation.\n",
    "\n",
    "It produces files \"flanksize_[flank size value].csv\" which provide a quantification of mutational concordance/homology between human and mouse as a function of flank size, which corresponds to the homology requirement of the codons flanking the site of the mutation.\n",
    "\n",
    "## Loading in required files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/samgould/anaconda3/lib/python3.7/site-packages/pandas/compat/_optional.py:138: UserWarning: Pandas requires version '2.7.0' or newer of 'numexpr' (version '2.6.8' currently installed).\n",
      "  warnings.warn(msg, UserWarning)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import csv\n",
    "from Bio import SeqIO\n",
    "import gzip\n",
    "from Bio.Seq import Seq\n",
    "import re\n",
    "import gffutils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns', 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/samgould/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3267: DtypeWarning: Columns (45,48,88) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Hugo_Symbol</th>\n",
       "      <th>Entrez_Gene_Id</th>\n",
       "      <th>Center</th>\n",
       "      <th>NCBI_Build</th>\n",
       "      <th>Chromosome</th>\n",
       "      <th>Start_Position</th>\n",
       "      <th>End_Position</th>\n",
       "      <th>Strand</th>\n",
       "      <th>Consequence</th>\n",
       "      <th>Variant_Classification</th>\n",
       "      <th>Variant_Type</th>\n",
       "      <th>Reference_Allele</th>\n",
       "      <th>Tumor_Seq_Allele1</th>\n",
       "      <th>Tumor_Seq_Allele2</th>\n",
       "      <th>dbSNP_RS</th>\n",
       "      <th>dbSNP_Val_Status</th>\n",
       "      <th>Tumor_Sample_Barcode</th>\n",
       "      <th>Matched_Norm_Sample_Barcode</th>\n",
       "      <th>Match_Norm_Seq_Allele1</th>\n",
       "      <th>Match_Norm_Seq_Allele2</th>\n",
       "      <th>Tumor_Validation_Allele1</th>\n",
       "      <th>Tumor_Validation_Allele2</th>\n",
       "      <th>Match_Norm_Validation_Allele1</th>\n",
       "      <th>Match_Norm_Validation_Allele2</th>\n",
       "      <th>Verification_Status</th>\n",
       "      <th>...</th>\n",
       "      <th>MOTIF_SCORE_CHANGE</th>\n",
       "      <th>PHENO</th>\n",
       "      <th>PICK</th>\n",
       "      <th>PUBMED</th>\n",
       "      <th>PolyPhen</th>\n",
       "      <th>SAS_MAF</th>\n",
       "      <th>SIFT</th>\n",
       "      <th>SOMATIC</th>\n",
       "      <th>SWISSPROT</th>\n",
       "      <th>SYMBOL</th>\n",
       "      <th>SYMBOL_SOURCE</th>\n",
       "      <th>TREMBL</th>\n",
       "      <th>TSL</th>\n",
       "      <th>Transcript</th>\n",
       "      <th>UNIPARC</th>\n",
       "      <th>VARIANT_CLASS</th>\n",
       "      <th>all_effects</th>\n",
       "      <th>amino_acid_change</th>\n",
       "      <th>cDNA_Change</th>\n",
       "      <th>cDNA_position</th>\n",
       "      <th>cdna_change</th>\n",
       "      <th>comments</th>\n",
       "      <th>n_depth</th>\n",
       "      <th>t_depth</th>\n",
       "      <th>transcript</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>BRCA2</td>\n",
       "      <td>675</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>GRCh37</td>\n",
       "      <td>13</td>\n",
       "      <td>32937315</td>\n",
       "      <td>32937315</td>\n",
       "      <td>+</td>\n",
       "      <td>splice_acceptor_variant</td>\n",
       "      <td>Splice_Site</td>\n",
       "      <td>SNP</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>rs81002874</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0029279-T01-IM6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>BRCA2</td>\n",
       "      <td>0</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>37</td>\n",
       "      <td>13</td>\n",
       "      <td>32914437</td>\n",
       "      <td>32914438</td>\n",
       "      <td>+</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>DEL</td>\n",
       "      <td>GT</td>\n",
       "      <td>GT</td>\n",
       "      <td>G</td>\n",
       "      <td>rs80359550</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0034227-T01-IM6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>MUTYH</td>\n",
       "      <td>4595</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>GRCh37</td>\n",
       "      <td>1</td>\n",
       "      <td>45798475</td>\n",
       "      <td>45798475</td>\n",
       "      <td>+</td>\n",
       "      <td>missense_variant</td>\n",
       "      <td>Missense_Mutation</td>\n",
       "      <td>SNP</td>\n",
       "      <td>T</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>rs34612342</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0030735-T01-IM6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>BRCA2</td>\n",
       "      <td>675</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>GRCh37</td>\n",
       "      <td>13</td>\n",
       "      <td>32893302</td>\n",
       "      <td>32893302</td>\n",
       "      <td>+</td>\n",
       "      <td>frameshift_variant</td>\n",
       "      <td>Frame_Shift_Ins</td>\n",
       "      <td>INS</td>\n",
       "      <td>T</td>\n",
       "      <td>T</td>\n",
       "      <td>GCCGGGCGCGGTGG</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0038798-T01-IM6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>BRCA1</td>\n",
       "      <td>0</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>37</td>\n",
       "      <td>17</td>\n",
       "      <td>41251824</td>\n",
       "      <td>41251825</td>\n",
       "      <td>+</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>DEL</td>\n",
       "      <td>TG</td>\n",
       "      <td>TG</td>\n",
       "      <td>T</td>\n",
       "      <td>rs80357872</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0030162-T01-IM6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>422817</th>\n",
       "      <td>SMARCA4</td>\n",
       "      <td>6597</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>GRCh37</td>\n",
       "      <td>19</td>\n",
       "      <td>11144132</td>\n",
       "      <td>11144132</td>\n",
       "      <td>+</td>\n",
       "      <td>missense_variant</td>\n",
       "      <td>Missense_Mutation</td>\n",
       "      <td>SNP</td>\n",
       "      <td>C</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0052864-T01-XS1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>422818</th>\n",
       "      <td>BRAF</td>\n",
       "      <td>673</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>GRCh37</td>\n",
       "      <td>7</td>\n",
       "      <td>140453149</td>\n",
       "      <td>140453149</td>\n",
       "      <td>+</td>\n",
       "      <td>missense_variant</td>\n",
       "      <td>Missense_Mutation</td>\n",
       "      <td>SNP</td>\n",
       "      <td>C</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>rs121913361</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0052867-T01-XS1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>422819</th>\n",
       "      <td>NRAS</td>\n",
       "      <td>4893</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>GRCh37</td>\n",
       "      <td>1</td>\n",
       "      <td>115258747</td>\n",
       "      <td>115258747</td>\n",
       "      <td>+</td>\n",
       "      <td>missense_variant</td>\n",
       "      <td>Missense_Mutation</td>\n",
       "      <td>SNP</td>\n",
       "      <td>C</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>rs121913237</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0052951-T01-XS1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>422820</th>\n",
       "      <td>TERT</td>\n",
       "      <td>7015</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>GRCh37</td>\n",
       "      <td>5</td>\n",
       "      <td>1295521</td>\n",
       "      <td>1295521</td>\n",
       "      <td>+</td>\n",
       "      <td>upstream_gene_variant</td>\n",
       "      <td>5'Flank</td>\n",
       "      <td>SNP</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>T</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0052951-T01-XS1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>422821</th>\n",
       "      <td>KRAS</td>\n",
       "      <td>3845</td>\n",
       "      <td>MSKCC</td>\n",
       "      <td>GRCh37</td>\n",
       "      <td>12</td>\n",
       "      <td>25398284</td>\n",
       "      <td>25398284</td>\n",
       "      <td>+</td>\n",
       "      <td>missense_variant</td>\n",
       "      <td>Missense_Mutation</td>\n",
       "      <td>SNP</td>\n",
       "      <td>C</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>P-0052952-T01-XS1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>422822 rows × 123 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       Hugo_Symbol  Entrez_Gene_Id Center NCBI_Build Chromosome  \\\n",
       "0            BRCA2             675  MSKCC     GRCh37         13   \n",
       "1            BRCA2               0  MSKCC         37         13   \n",
       "2            MUTYH            4595  MSKCC     GRCh37          1   \n",
       "3            BRCA2             675  MSKCC     GRCh37         13   \n",
       "4            BRCA1               0  MSKCC         37         17   \n",
       "...            ...             ...    ...        ...        ...   \n",
       "422817     SMARCA4            6597  MSKCC     GRCh37         19   \n",
       "422818        BRAF             673  MSKCC     GRCh37          7   \n",
       "422819        NRAS            4893  MSKCC     GRCh37          1   \n",
       "422820        TERT            7015  MSKCC     GRCh37          5   \n",
       "422821        KRAS            3845  MSKCC     GRCh37         12   \n",
       "\n",
       "        Start_Position  End_Position Strand              Consequence  \\\n",
       "0             32937315      32937315      +  splice_acceptor_variant   \n",
       "1             32914437      32914438      +                      NaN   \n",
       "2             45798475      45798475      +         missense_variant   \n",
       "3             32893302      32893302      +       frameshift_variant   \n",
       "4             41251824      41251825      +                      NaN   \n",
       "...                ...           ...    ...                      ...   \n",
       "422817        11144132      11144132      +         missense_variant   \n",
       "422818       140453149     140453149      +         missense_variant   \n",
       "422819       115258747     115258747      +         missense_variant   \n",
       "422820         1295521       1295521      +    upstream_gene_variant   \n",
       "422821        25398284      25398284      +         missense_variant   \n",
       "\n",
       "       Variant_Classification Variant_Type Reference_Allele Tumor_Seq_Allele1  \\\n",
       "0                 Splice_Site          SNP                G                 G   \n",
       "1                         NaN          DEL               GT                GT   \n",
       "2           Missense_Mutation          SNP                T                 T   \n",
       "3             Frame_Shift_Ins          INS                T                 T   \n",
       "4                         NaN          DEL               TG                TG   \n",
       "...                       ...          ...              ...               ...   \n",
       "422817      Missense_Mutation          SNP                C                 C   \n",
       "422818      Missense_Mutation          SNP                C                 C   \n",
       "422819      Missense_Mutation          SNP                C                 C   \n",
       "422820                5'Flank          SNP                A                 A   \n",
       "422821      Missense_Mutation          SNP                C                 C   \n",
       "\n",
       "       Tumor_Seq_Allele2     dbSNP_RS  dbSNP_Val_Status Tumor_Sample_Barcode  \\\n",
       "0                      C   rs81002874               NaN    P-0029279-T01-IM6   \n",
       "1                      G   rs80359550               NaN    P-0034227-T01-IM6   \n",
       "2                      C   rs34612342               NaN    P-0030735-T01-IM6   \n",
       "3         GCCGGGCGCGGTGG          NaN               NaN    P-0038798-T01-IM6   \n",
       "4                      T   rs80357872               NaN    P-0030162-T01-IM6   \n",
       "...                  ...          ...               ...                  ...   \n",
       "422817                 G          NaN               NaN    P-0052864-T01-XS1   \n",
       "422818                 G  rs121913361               NaN    P-0052867-T01-XS1   \n",
       "422819                 T  rs121913237               NaN    P-0052951-T01-XS1   \n",
       "422820                 T          NaN               NaN    P-0052951-T01-XS1   \n",
       "422821                 A          NaN               NaN    P-0052952-T01-XS1   \n",
       "\n",
       "        Matched_Norm_Sample_Barcode  Match_Norm_Seq_Allele1  \\\n",
       "0                               NaN                     NaN   \n",
       "1                               NaN                     NaN   \n",
       "2                               NaN                     NaN   \n",
       "3                               NaN                     NaN   \n",
       "4                               NaN                     NaN   \n",
       "...                             ...                     ...   \n",
       "422817                          NaN                     NaN   \n",
       "422818                          NaN                     NaN   \n",
       "422819                          NaN                     NaN   \n",
       "422820                          NaN                     NaN   \n",
       "422821                          NaN                     NaN   \n",
       "\n",
       "        Match_Norm_Seq_Allele2  Tumor_Validation_Allele1  \\\n",
       "0                          NaN                       NaN   \n",
       "1                          NaN                       NaN   \n",
       "2                          NaN                       NaN   \n",
       "3                          NaN                       NaN   \n",
       "4                          NaN                       NaN   \n",
       "...                        ...                       ...   \n",
       "422817                     NaN                       NaN   \n",
       "422818                     NaN                       NaN   \n",
       "422819                     NaN                       NaN   \n",
       "422820                     NaN                       NaN   \n",
       "422821                     NaN                       NaN   \n",
       "\n",
       "        Tumor_Validation_Allele2  Match_Norm_Validation_Allele1  \\\n",
       "0                            NaN                            NaN   \n",
       "1                            NaN                            NaN   \n",
       "2                            NaN                            NaN   \n",
       "3                            NaN                            NaN   \n",
       "4                            NaN                            NaN   \n",
       "...                          ...                            ...   \n",
       "422817                       NaN                            NaN   \n",
       "422818                       NaN                            NaN   \n",
       "422819                       NaN                            NaN   \n",
       "422820                       NaN                            NaN   \n",
       "422821                       NaN                            NaN   \n",
       "\n",
       "        Match_Norm_Validation_Allele2  Verification_Status  ...  \\\n",
       "0                                 NaN                  NaN  ...   \n",
       "1                                 NaN                  NaN  ...   \n",
       "2                                 NaN                  NaN  ...   \n",
       "3                                 NaN                  NaN  ...   \n",
       "4                                 NaN                  NaN  ...   \n",
       "...                               ...                  ...  ...   \n",
       "422817                            NaN                  NaN  ...   \n",
       "422818                            NaN                  NaN  ...   \n",
       "422819                            NaN                  NaN  ...   \n",
       "422820                            NaN                  NaN  ...   \n",
       "422821                            NaN                  NaN  ...   \n",
       "\n",
       "       MOTIF_SCORE_CHANGE PHENO  PICK  PUBMED  PolyPhen SAS_MAF  SIFT  \\\n",
       "0                     NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "1                     NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "2                     NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "3                     NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "4                     NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "...                   ...   ...   ...     ...       ...     ...   ...   \n",
       "422817                NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "422818                NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "422819                NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "422820                NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "422821                NaN   NaN   NaN     NaN       NaN     NaN   NaN   \n",
       "\n",
       "        SOMATIC  SWISSPROT  SYMBOL  SYMBOL_SOURCE  TREMBL TSL Transcript  \\\n",
       "0           NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "1           NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "2           NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "3           NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "4           NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "...         ...        ...     ...            ...     ...  ..        ...   \n",
       "422817      NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "422818      NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "422819      NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "422820      NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "422821      NaN        NaN     NaN            NaN     NaN NaN        NaN   \n",
       "\n",
       "       UNIPARC VARIANT_CLASS all_effects  amino_acid_change cDNA_Change  \\\n",
       "0          NaN           NaN         NaN                NaN         NaN   \n",
       "1          NaN           NaN         NaN                NaN         NaN   \n",
       "2          NaN           NaN         NaN                NaN         NaN   \n",
       "3          NaN           NaN         NaN                NaN         NaN   \n",
       "4          NaN           NaN         NaN                NaN         NaN   \n",
       "...        ...           ...         ...                ...         ...   \n",
       "422817     NaN           NaN         NaN                NaN         NaN   \n",
       "422818     NaN           NaN         NaN                NaN         NaN   \n",
       "422819     NaN           NaN         NaN                NaN         NaN   \n",
       "422820     NaN           NaN         NaN                NaN         NaN   \n",
       "422821     NaN           NaN         NaN                NaN         NaN   \n",
       "\n",
       "        cDNA_position cdna_change  comments  n_depth t_depth  transcript  \n",
       "0                 NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "1                 NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "2                 NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "3                 NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "4                 NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "...               ...         ...       ...      ...     ...         ...  \n",
       "422817            NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "422818            NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "422819            NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "422820            NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "422821            NaN         NaN       NaN      NaN     NaN         NaN  \n",
       "\n",
       "[422822 rows x 123 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filepath = '/Volumes/Sam_G_SSD/2020-06-16-MSK-IMPACT_EDITED.txt'\n",
    "impact_data = pd.read_csv(filepath, sep='\\t')\n",
    "impact_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>gene</th>\n",
       "      <th>human gene name</th>\n",
       "      <th>mouse gene name</th>\n",
       "      <th>mouse id</th>\n",
       "      <th>mouse id version</th>\n",
       "      <th>mouse transcript</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ABL1</td>\n",
       "      <td>ABL1</td>\n",
       "      <td>Abl1</td>\n",
       "      <td>ENSMUSG00000026842</td>\n",
       "      <td>ENSMUSG00000026842.16</td>\n",
       "      <td>ENSMUST00000028190.12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AC004906.3</td>\n",
       "      <td>NONE</td>\n",
       "      <td>NONE</td>\n",
       "      <td>NONE</td>\n",
       "      <td>NONE</td>\n",
       "      <td>NONE</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AC008738.1</td>\n",
       "      <td>NONE</td>\n",
       "      <td>NONE</td>\n",
       "      <td>NONE</td>\n",
       "      <td>NONE</td>\n",
       "      <td>NONE</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ACTG1</td>\n",
       "      <td>ACTG1</td>\n",
       "      <td>Actg1</td>\n",
       "      <td>ENSMUSG00000062825</td>\n",
       "      <td>ENSMUSG00000062825.15</td>\n",
       "      <td>ENSMUST00000071555.12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ACVR1</td>\n",
       "      <td>ACVR1</td>\n",
       "      <td>Acvr1</td>\n",
       "      <td>ENSMUSG00000026836</td>\n",
       "      <td>ENSMUSG00000026836.15</td>\n",
       "      <td>ENSMUST00000056376.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>589</th>\n",
       "      <td>XRCC2</td>\n",
       "      <td>XRCC2</td>\n",
       "      <td>Xrcc2</td>\n",
       "      <td>ENSMUSG00000028933</td>\n",
       "      <td>ENSMUSG00000028933.11</td>\n",
       "      <td>ENSMUST00000030773.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>590</th>\n",
       "      <td>YAP1</td>\n",
       "      <td>YAP1</td>\n",
       "      <td>Yap1</td>\n",
       "      <td>ENSMUSG00000053110</td>\n",
       "      <td>ENSMUSG00000053110.13</td>\n",
       "      <td>ENSMUST00000086580.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>591</th>\n",
       "      <td>YES1</td>\n",
       "      <td>YES1</td>\n",
       "      <td>Yes1</td>\n",
       "      <td>ENSMUSG00000014932</td>\n",
       "      <td>ENSMUSG00000014932.15</td>\n",
       "      <td>ENSMUST00000168707.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>592</th>\n",
       "      <td>ZFHX3</td>\n",
       "      <td>ZFHX3</td>\n",
       "      <td>Zfhx3</td>\n",
       "      <td>ENSMUSG00000038872</td>\n",
       "      <td>ENSMUSG00000038872.10</td>\n",
       "      <td>ENSMUST00000043896.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>593</th>\n",
       "      <td>ZRSR2</td>\n",
       "      <td>ZRSR2</td>\n",
       "      <td>Zrsr1</td>\n",
       "      <td>ENSMUSG00000044068</td>\n",
       "      <td>ENSMUSG00000044068.7</td>\n",
       "      <td>ENSMUST00000049506.7</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>594 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           gene human gene name mouse gene name            mouse id  \\\n",
       "0          ABL1            ABL1            Abl1  ENSMUSG00000026842   \n",
       "1    AC004906.3            NONE            NONE                NONE   \n",
       "2    AC008738.1            NONE            NONE                NONE   \n",
       "3         ACTG1           ACTG1           Actg1  ENSMUSG00000062825   \n",
       "4         ACVR1           ACVR1           Acvr1  ENSMUSG00000026836   \n",
       "..          ...             ...             ...                 ...   \n",
       "589       XRCC2           XRCC2           Xrcc2  ENSMUSG00000028933   \n",
       "590        YAP1            YAP1            Yap1  ENSMUSG00000053110   \n",
       "591        YES1            YES1            Yes1  ENSMUSG00000014932   \n",
       "592       ZFHX3           ZFHX3           Zfhx3  ENSMUSG00000038872   \n",
       "593       ZRSR2           ZRSR2           Zrsr1  ENSMUSG00000044068   \n",
       "\n",
       "          mouse id version       mouse transcript  \n",
       "0    ENSMUSG00000026842.16  ENSMUST00000028190.12  \n",
       "1                     NONE                   NONE  \n",
       "2                     NONE                   NONE  \n",
       "3    ENSMUSG00000062825.15  ENSMUST00000071555.12  \n",
       "4    ENSMUSG00000026836.15  ENSMUST00000056376.11  \n",
       "..                     ...                    ...  \n",
       "589  ENSMUSG00000028933.11  ENSMUST00000030773.11  \n",
       "590  ENSMUSG00000053110.13  ENSMUST00000086580.11  \n",
       "591  ENSMUSG00000014932.15   ENSMUST00000168707.5  \n",
       "592  ENSMUSG00000038872.10   ENSMUST00000043896.9  \n",
       "593   ENSMUSG00000044068.7   ENSMUST00000049506.7  \n",
       "\n",
       "[594 rows x 6 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#homology dataframe\n",
    "\n",
    "homology_df = np.load('/Volumes/Sam_G_SSD/homology_table.npy', allow_pickle=True)\n",
    "homology_df = pd.DataFrame(homology_df,columns=['gene','human gene name','mouse gene name','mouse id','mouse id version','mouse transcript'])\n",
    "\n",
    "homology_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>gene</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>transcript_id</th>\n",
       "      <th>chrom</th>\n",
       "      <th>gene_start</th>\n",
       "      <th>gene_end</th>\n",
       "      <th>transcript_start</th>\n",
       "      <th>transcript_end</th>\n",
       "      <th>strand</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ABL1</td>\n",
       "      <td>ENSG00000097007.13</td>\n",
       "      <td>ENST00000318560.5</td>\n",
       "      <td>chr9</td>\n",
       "      <td>133589333</td>\n",
       "      <td>133763062</td>\n",
       "      <td>133710453</td>\n",
       "      <td>133763062</td>\n",
       "      <td>+</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AC004906.3</td>\n",
       "      <td>ENSG00000237286.1</td>\n",
       "      <td>ENST00000423194.1</td>\n",
       "      <td>chr7</td>\n",
       "      <td>2983669</td>\n",
       "      <td>2986725</td>\n",
       "      <td>2983669</td>\n",
       "      <td>2986725</td>\n",
       "      <td>+</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AC008738.1</td>\n",
       "      <td>ENSG00000230259.2</td>\n",
       "      <td>ENST00000425420.2</td>\n",
       "      <td>chr19</td>\n",
       "      <td>33790853</td>\n",
       "      <td>33793430</td>\n",
       "      <td>33790853</td>\n",
       "      <td>33793430</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ACTG1</td>\n",
       "      <td>ENSG00000184009.5</td>\n",
       "      <td>ENST00000575842.1</td>\n",
       "      <td>chr17</td>\n",
       "      <td>79476997</td>\n",
       "      <td>79490873</td>\n",
       "      <td>79477015</td>\n",
       "      <td>79479807</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ACVR1</td>\n",
       "      <td>ENSG00000115170.9</td>\n",
       "      <td>ENST00000263640.3</td>\n",
       "      <td>chr2</td>\n",
       "      <td>158592958</td>\n",
       "      <td>158732374</td>\n",
       "      <td>158592958</td>\n",
       "      <td>158731623</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>589</th>\n",
       "      <td>XRCC2</td>\n",
       "      <td>ENSG00000196584.2</td>\n",
       "      <td>ENST00000359321.1</td>\n",
       "      <td>chr7</td>\n",
       "      <td>152341864</td>\n",
       "      <td>152373250</td>\n",
       "      <td>152343589</td>\n",
       "      <td>152373250</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>590</th>\n",
       "      <td>YAP1</td>\n",
       "      <td>ENSG00000137693.9</td>\n",
       "      <td>ENST00000282441.5</td>\n",
       "      <td>chr11</td>\n",
       "      <td>101981192</td>\n",
       "      <td>102104154</td>\n",
       "      <td>101981192</td>\n",
       "      <td>102104154</td>\n",
       "      <td>+</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>591</th>\n",
       "      <td>YES1</td>\n",
       "      <td>ENSG00000176105.9</td>\n",
       "      <td>ENST00000314574.4</td>\n",
       "      <td>chr18</td>\n",
       "      <td>721588</td>\n",
       "      <td>812547</td>\n",
       "      <td>721748</td>\n",
       "      <td>812239</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>592</th>\n",
       "      <td>ZFHX3</td>\n",
       "      <td>ENSG00000140836.10</td>\n",
       "      <td>ENST00000268489.5</td>\n",
       "      <td>chr16</td>\n",
       "      <td>72816784</td>\n",
       "      <td>73093597</td>\n",
       "      <td>72816784</td>\n",
       "      <td>73082274</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>593</th>\n",
       "      <td>ZRSR2</td>\n",
       "      <td>ENSG00000169249.8</td>\n",
       "      <td>ENST00000307771.7</td>\n",
       "      <td>chrX</td>\n",
       "      <td>15808595</td>\n",
       "      <td>15841383</td>\n",
       "      <td>15808595</td>\n",
       "      <td>15841383</td>\n",
       "      <td>+</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>594 rows × 9 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           gene             gene_id      transcript_id  chrom  gene_start  \\\n",
       "0          ABL1  ENSG00000097007.13  ENST00000318560.5   chr9   133589333   \n",
       "1    AC004906.3   ENSG00000237286.1  ENST00000423194.1   chr7     2983669   \n",
       "2    AC008738.1   ENSG00000230259.2  ENST00000425420.2  chr19    33790853   \n",
       "3         ACTG1   ENSG00000184009.5  ENST00000575842.1  chr17    79476997   \n",
       "4         ACVR1   ENSG00000115170.9  ENST00000263640.3   chr2   158592958   \n",
       "..          ...                 ...                ...    ...         ...   \n",
       "589       XRCC2   ENSG00000196584.2  ENST00000359321.1   chr7   152341864   \n",
       "590        YAP1   ENSG00000137693.9  ENST00000282441.5  chr11   101981192   \n",
       "591        YES1   ENSG00000176105.9  ENST00000314574.4  chr18      721588   \n",
       "592       ZFHX3  ENSG00000140836.10  ENST00000268489.5  chr16    72816784   \n",
       "593       ZRSR2   ENSG00000169249.8  ENST00000307771.7   chrX    15808595   \n",
       "\n",
       "      gene_end  transcript_start  transcript_end strand  \n",
       "0    133763062         133710453       133763062      +  \n",
       "1      2986725           2983669         2986725      +  \n",
       "2     33793430          33790853        33793430      -  \n",
       "3     79490873          79477015        79479807      -  \n",
       "4    158732374         158592958       158731623      -  \n",
       "..         ...               ...             ...    ...  \n",
       "589  152373250         152343589       152373250      -  \n",
       "590  102104154         101981192       102104154      +  \n",
       "591     812547            721748          812239      -  \n",
       "592   73093597          72816784        73082274      -  \n",
       "593   15841383          15808595        15841383      +  \n",
       "\n",
       "[594 rows x 9 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filename1 = '/Users/samgould/Desktop/FSR Lab/2022-03-17/gene_info.csv'\n",
    "df1 = pd.read_csv(filename1)\n",
    "df1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#loading in annotation databases for human and mouse\n",
    "file = '/Volumes/Sam_G_SSD/gencode_v19.db'\n",
    "db = gffutils.FeatureDB(file)\n",
    "\n",
    "file_mouse = '/Volumes/Sam_G_SSD/GRCm38.p6 (mouse)/gencode_vM25.db'\n",
    "db_mouse = gffutils.FeatureDB(file_mouse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "#loading in necessary genes for human\n",
    "path = '/Volumes/Sam_G_SSD/human genome GrCh37 IMPACT genes/'\n",
    "impact_genes = np.load(path + 'human_impact_genes_plusminus5000.npy', allow_pickle=True)\n",
    "unique_genes = np.load(path + 'human_impact_genes_NAMES.npy', allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#loading in mouse genes\n",
    "#loading in necessary genes for human\n",
    "path = '/Volumes/Sam_G_SSD/mouse genome GRCm38.p6 IMPACT genes/'\n",
    "mouse_genes = np.load(path + 'mouse_impact_genes_plusminus5000.npy', allow_pickle=True)\n",
    "mouse_gene_names = np.load(path + 'mouse_impact_genes_NAMES.npy', allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'+'"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "db_mouse['ENSMUSG00000026842.16'].strand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#loading in orthologous protein alignmnet\n",
    "\n",
    "path = '/Volumes/Sam_G_SSD/human_mouse_alignments/'\n",
    "human_prot_align = np.load(path+'human_alignment_idx.npy', allow_pickle=True)\n",
    "mouse_prot_align = np.load(path+'mouse_alignment_idx.npy', allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#loading in mapping between codons and amino acids\n",
    "path = '/Volumes/Sam_G_SSD/human_mouse_alignments/'\n",
    "\n",
    "human_codon_locations = np.load(path+'human_codon_locations.npy', allow_pickle=True)\n",
    "human_codon_seqs = np.load(path+'human_codon_seqs.npy', allow_pickle=True)\n",
    "human_aa = np.load(path+'human_aa.npy', allow_pickle=True)\n",
    "\n",
    "mouse_codon_locations = np.load(path+'mouse_codon_locations.npy', allow_pickle=True)\n",
    "mouse_codon_seqs = np.load(path+'mouse_codon_seqs.npy', allow_pickle=True)\n",
    "mouse_aa = np.load(path+'mouse_aa.npy', allow_pickle=True)\n",
    "\n",
    "excluded_genes = []\n",
    "for i in range(len(mouse_aa)):\n",
    "    \n",
    "    if len(mouse_aa[i])==1:\n",
    "        excluded_genes.append(unique_genes[i])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Putting everything together to determine if homologous mutation can be modeled\n",
    "\n",
    "- In the below script, I am able to determine whether a given mutation falls in a region of alignment (i.e. homology) at varying stringencies of homologous flanking region.\n",
    "- However, the script is currently incapable of modeling the effects of each mutation in human/mouse. I do, however, record if both the amino acid and DNA sequence are conserved, but this does not take into account potentially synonymous mutations with differing DNA sequences and conserved amino acid sequences. In short, it needs to be fixed, but it is not the current focus of this analysis.\n",
    "    - Note: the main issue preventing it from being fixed easily is that all mutations are reported on the + strand, but not all genes are transcribed in the + direction (and this directionality doesn't need to match between species). Essentially, it is a difficult indexing problem that could be fixed if enough time were devoted to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "code_folding": [
     23,
     137
    ]
   },
   "outputs": [],
   "source": [
    "#function that takes in original codon(s) and mutation types, and spits out resulting codon sequence\n",
    "#see above for errors associated with it\n",
    "\n",
    "def mutation_modeling(h_codon, m_codon, codon_locs, codon_idx,human_aa_idx, mouse_aa_idx, within_codon_idx, mutation, gene_idx):\n",
    "    variant_type = mutation['Variant_Type'].values[0]\n",
    "    start = mutation['Start_Position'].values[0]\n",
    "    end = mutation['End_Position'].values[0]\n",
    "    \n",
    "    ref_allele = mutation['Reference_Allele'].values[0]\n",
    "    mut_allele = mutation['Tumor_Seq_Allele2'].values[0]\n",
    "    \n",
    "    #for recording \n",
    "    aa_concordant = 0\n",
    "    dna_concordant = 0\n",
    "    ref_aa = human_aa[gene_idx][codon_idx]\n",
    "    mut_aa_h = 'none'\n",
    "    mut_aa_m = 'none'\n",
    "    \n",
    "    #defining which strand the gene falls on in mouse and human\n",
    "    tx = homology_df[homology_df['gene']==gene_name]['mouse transcript'].values[0]\n",
    "    strand_m = db_mouse[tx].strand\n",
    "    strand_h = df1[df1['gene']==gene_name]['strand'].values[0]\n",
    "    \n",
    "    #modeling mutations\n",
    "    \n",
    "    if variant_type=='DEL':\n",
    "        \n",
    "        end_ind = np.where(np.array(codon_locs) == end) #where does the mutation fall (which codon)\n",
    "        codon_idx_end = end_ind[0][0] #which codon is it\n",
    "        within_codon_idx_end = end_ind[1][0]\n",
    "        \n",
    "        diff = abs(codon_idx_end - codon_idx)\n",
    "        \n",
    "        seq_ref_h = sum(human_codon_seqs[gene_idx][human_aa_idx:human_aa_idx+diff+1], Seq(''))\n",
    "        #need to be mindful of strand...\n",
    "        aa_ref_h = sum(human_aa[gene_idx][human_aa_idx:human_aa_idx+diff+1], Seq(''))\n",
    "        \n",
    "        ref_aa = aa_ref_h #record reference amino acids that are spanned\n",
    "        \n",
    "        aa_ref_m = sum(mouse_aa[gene_idx][mouse_aa_idx:mouse_aa_idx+diff+1], Seq(''))\n",
    "        seq_ref_m = sum(mouse_codon_seqs[gene_idx][mouse_aa_idx:mouse_aa_idx+diff+1], Seq(''))\n",
    "        \n",
    "        tx = homology_df[homology_df['gene']==gene_name]['mouse transcript'].values[0]\n",
    "        strand_m = db_mouse[tx].strand\n",
    "        if strand_m=='-': #if minus strand take complement\n",
    "            seq_ref_m = seq_ref_m.complement()\n",
    "        elif strand_m=='+':\n",
    "            seq_ref_m = seq_ref_m\n",
    " \n",
    "        #if aa level is concordant records\n",
    "        if aa_ref_h == aa_ref_m:\n",
    "            aa_concordant = 1\n",
    "        if seq_ref_h == seq_ref_m:\n",
    "            dna_concordant = 1\n",
    "    \n",
    "      \n",
    "    elif variant_type=='INS':\n",
    "        \n",
    "        end_ind = np.where(np.array(codon_locs) == end) #where does the mutation fall (which codon)\n",
    "        codon_idx_end = end_ind[0][0] #which codon is it\n",
    "        within_codon_idx_end = end_ind[1][0]\n",
    "\n",
    "        diff = abs(codon_idx_end - codon_idx)\n",
    "\n",
    "        tx = homology_df[homology_df['gene']==gene_name]['mouse transcript'].values[0]\n",
    "        strand_m = db_mouse[tx].strand\n",
    "        strand_h = df1[df1['gene']==gene_name]['strand'].values[0]    \n",
    "\n",
    "        seq_ref_h = sum(human_codon_seqs[gene_idx][human_aa_idx:human_aa_idx+diff+1], Seq(''))\n",
    "                #need to be mindful of strand...        \n",
    "        aa_ref_h = sum(human_aa[gene_idx][human_aa_idx:human_aa_idx+diff+1], Seq(''))\n",
    "\n",
    "        seq_ref_m = sum(mouse_codon_seqs[gene_idx][mouse_aa_idx:mouse_aa_idx+diff+1], Seq(''))\n",
    "\n",
    "        within_codon_idx_h = within_codon_idx\n",
    "        within_codon_idx_end_h = within_codon_idx_end\n",
    "\n",
    "\n",
    "        if strand_h=='-':\n",
    "            within_codon_idx_h = 2- within_codon_idx\n",
    "            within_codon_idx_end_h = 2 - within_codon_idx_end\n",
    "\n",
    "        #ref_aa = aa_ref_h #record reference amino acids that are spanned\n",
    "        ind_2 = np.where(np.array(codon_locs) == codon_idx+diff)\n",
    "        if len(ind_2[0])>0:\n",
    "            if strand_h==strand_m:\n",
    "                mutant_seq_h = h_codon[0:within_codon_idx_h] + mut_allele + human_codon_seqs[gene_idx][codon_idx+diff][within_codon_idx_h:]\n",
    "                mutant_seq_m = m_codon[0:within_codon_idx_h] + mut_allele + mouse_codon_seqs[gene_idx][mouse_aa_idx+diff][within_codon_idx_h:]\n",
    "\n",
    "            else:\n",
    "                mutant_seq_h = h_codon[0:within_codon_idx_h] + mut_allele + human_codon_seqs[gene_idx][codon_idx+diff][within_codon_idx_h:]\n",
    "                mutant_seq_m = m_codon[0:within_codon_idx_h] + Seq(mut_allele).complement() + mouse_codon_seqs[gene_idx][mouse_aa_idx+diff][within_codon_idx_h:]\n",
    "\n",
    "\n",
    "            right_flank_idx_h = min(human_aa_idx+1, max(human_prot_align[gene_idx]))\n",
    "            right_flank_idx_m = min(mouse_aa_idx+1, max(mouse_prot_align[gene_idx]))\n",
    "            mutant_seq_h_flanked = human_codon_seqs[gene_idx][human_aa_idx-1] + mutant_seq_h + human_codon_seqs[gene_idx][right_flank_idx_h]\n",
    "            mutant_seq_m_flanked = mouse_codon_seqs[gene_idx][mouse_aa_idx-1] + mutant_seq_m + mouse_codon_seqs[gene_idx][right_flank_idx_m]\n",
    "\n",
    "            strand_h = df1[df1['gene']==gene_name]['strand'].values[0]\n",
    "            if strand_h=='-': #if minus strand\n",
    "                mutant_aa_h = mutant_seq_h_flanked.complement().transcribe().translate()\n",
    "                ref_aa = seq_ref_h.complement().transcribe().translate()\n",
    "            elif strand_h=='+':\n",
    "                mutant_aa_h = mutant_seq_h_flanked.transcribe().translate()\n",
    "                ref_aa = seq_ref_h.transcribe().translate()\n",
    "\n",
    "            tx = homology_df[homology_df['gene']==gene_name]['mouse transcript'].values[0]\n",
    "            strand_m = db_mouse[tx].strand\n",
    "            if strand_m=='-': #if minus strand take complement\n",
    "                mutant_aa_m = mutant_seq_m_flanked.complement().transcribe().translate()\n",
    "            elif strand_m=='+':\n",
    "                mutant_aa_m = mutant_seq_m_flanked.transcribe().translate()\n",
    "\n",
    "            mut_aa_h = mutant_aa_h\n",
    "            mut_aa_m = mutant_aa_m\n",
    "\n",
    "            if mut_aa_h==mut_aa_m:\n",
    "                aa_concordant=1\n",
    "            else:\n",
    "                aa_concordant=0        \n",
    "\n",
    "            if mutant_aa_h==mutant_aa_m:\n",
    "                aa_concordant=1\n",
    "            else:\n",
    "                aa_concordant=0\n",
    "\n",
    "            #need to double check that this is correct...\n",
    "            if strand_h==strand_m:\n",
    "                if mutant_seq_m_flanked == mutant_seq_h_flanked:\n",
    "                    dna_concordant=1\n",
    "                else:\n",
    "                    dna_concordant=0\n",
    "\n",
    "            else:\n",
    "                if mutant_seq_m_flanked.complement() == mutant_seq_h_flanked:\n",
    "                    dna_concordant=1\n",
    "                else:\n",
    "                    dna_concordant=0\n",
    "        \n",
    "        \n",
    "        #else:continue\n",
    "            \n",
    "            \n",
    "    else: #SNPs, ONPs, DNPs\n",
    "        end_ind = np.where(np.array(codon_locs) == end) #where does the mutation fall (which codon)\n",
    "        codon_idx_end = end_ind[0][0] #which codon is it\n",
    "        within_codon_idx_end = end_ind[1][0]\n",
    "\n",
    "        diff = abs(codon_idx_end - codon_idx)\n",
    "\n",
    "        tx = homology_df[homology_df['gene']==gene_name]['mouse transcript'].values[0]\n",
    "        strand_m = db_mouse[tx].strand\n",
    "        strand_h = df1[df1['gene']==gene_name]['strand'].values[0]    \n",
    "\n",
    "        seq_ref_h = sum(human_codon_seqs[gene_idx][human_aa_idx:human_aa_idx+diff+1], Seq(''))\n",
    "                #need to be mindful of strand...        \n",
    "        aa_ref_h = sum(human_aa[gene_idx][human_aa_idx:human_aa_idx+diff+1], Seq(''))\n",
    "\n",
    "        seq_ref_m = sum(mouse_codon_seqs[gene_idx][mouse_aa_idx:mouse_aa_idx+diff+1], Seq(''))\n",
    "\n",
    "        within_codon_idx_h = within_codon_idx\n",
    "        within_codon_idx_end_h = within_codon_idx_end\n",
    "\n",
    "\n",
    "        if strand_h=='-': #flip indexing if reverse strand\n",
    "            within_codon_idx_h = 2- within_codon_idx\n",
    "            within_codon_idx_end_h = 2 - within_codon_idx_end\n",
    "\n",
    "        #ref_aa = aa_ref_h #record reference amino acids that are spanned\n",
    "        if strand_h==strand_m:\n",
    "            mutant_seq_h = h_codon[0:within_codon_idx_h] + mut_allele + human_codon_seqs[gene_idx][codon_idx+diff][within_codon_idx_end_h+1:]\n",
    "            mutant_seq_m = m_codon[0:within_codon_idx_h] + mut_allele + mouse_codon_seqs[gene_idx][mouse_aa_idx+diff][within_codon_idx_end_h+1:]\n",
    "\n",
    "        else:\n",
    "            mutant_seq_h = h_codon[0:within_codon_idx_h] + mut_allele + human_codon_seqs[gene_idx][codon_idx+diff][within_codon_idx_end_h+1:]\n",
    "            mutant_seq_m = m_codon[0:within_codon_idx_h] + Seq(mut_allele).complement() + mouse_codon_seqs[gene_idx][mouse_aa_idx+diff][within_codon_idx_end_h+1:]\n",
    "\n",
    "\n",
    "        #need to double check that this is correct...\n",
    "\n",
    "        if strand_h=='-':\n",
    "            mutant_aa_h = mutant_seq_h.complement().transcribe().translate()\n",
    "            ref_aa = seq_ref_h.complement().transcribe().translate()\n",
    "        elif strand_h=='+':\n",
    "            mutant_aa_h = mutant_seq_h.transcribe().translate()\n",
    "            ref_aa = seq_ref_h.transcribe().translate()\n",
    "            \n",
    "        if strand_m=='-':\n",
    "            mutant_aa_m = mutant_seq_m.complement().transcribe().translate()\n",
    "        elif strand_m=='+':\n",
    "            mutant_aa_m = mutant_seq_m.transcribe().translate()\n",
    "\n",
    "        if strand_h==strand_m:\n",
    "            if mutant_seq_m == mutant_seq_h:\n",
    "                dna_concordant=1\n",
    "            else:\n",
    "                dna_concordant=0\n",
    "\n",
    "        else:\n",
    "            if mutant_seq_m.complement() == mutant_seq_h:\n",
    "                dna_concordant=1\n",
    "            else:\n",
    "                dna_concordant=0\n",
    "\n",
    "        if mutant_aa_h==mutant_aa_m:\n",
    "            aa_concordant=1\n",
    "        else:\n",
    "            aa_concordant=0\n",
    "\n",
    "    \n",
    "        mut_aa_h = mutant_aa_h\n",
    "        mut_aa_m = mutant_aa_m\n",
    "    \n",
    "    #all of the amino acids at the target site are concordanta t this point\n",
    "    #this is referring to the CONSEQUENCE...\n",
    "    return aa_concordant, dna_concordant, ref_aa, mut_aa_h, mut_aa_m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "#flank_size=2\n",
    "def ortholog_PE(flank_size):\n",
    "    unique_gene_list = list(unique_genes)\n",
    "    num_mutations =  len(impact_data)\n",
    "    #num_mutations=10000\n",
    "    excluded_mutations = np.zeros(num_mutations)\n",
    "    non_coding_mutations = np.zeros(num_mutations)\n",
    "\n",
    "    coding_mutations = np.zeros(num_mutations)\n",
    "    homologous = np.zeros(num_mutations) #refers to homologous at site of mutation (notconsequenece)\n",
    "    non_homologous = np.zeros(num_mutations)\n",
    "\n",
    "    #recording mutational concsequences\n",
    "    aa_concordant_record = np.zeros(num_mutations)\n",
    "    dna_concordant_record = np.zeros(num_mutations)\n",
    "\n",
    "    #only recording these for homologous mutations\n",
    "    original_aa = []\n",
    "    new_aa_h = []\n",
    "    new_aa_m = []\n",
    "\n",
    "    #iterate over the mutations\n",
    "    for i in range(num_mutations):\n",
    "        \n",
    "        mutation = impact_data.iloc[[i]]\n",
    "        gene_name = mutation['Hugo_Symbol'].values[0]\n",
    "        gene_idx = unique_gene_list.index(gene_name)\n",
    "\n",
    "        start = mutation['Start_Position'].values[0]\n",
    "        end = mutation['End_Position'].values[0]\n",
    "\n",
    "        codon_locs = human_codon_locations[gene_idx]\n",
    "        ind = np.where(np.array(codon_locs) == start) #where does the mutation fall (which codon)\n",
    "        end_ind = np.where(np.array(codon_locs) == end) #where does the mutation fall (which codon)\n",
    "\n",
    "\n",
    "        #first check if the mutation falls in an excluded gene\n",
    "        if gene_name in excluded_genes:\n",
    "            excluded_mutations[i]=1\n",
    "            original_aa.append('none')\n",
    "            new_aa_h.append('none')\n",
    "            new_aa_m.append('none')\n",
    "\n",
    "        #then check if it falls in a coding sequence\n",
    "        elif len(ind[0])==0:\n",
    "            non_coding_mutations[i]=1\n",
    "            original_aa.append('none')\n",
    "            new_aa_h.append('none')\n",
    "            new_aa_m.append('none')\n",
    "\n",
    "            #then check if it falls in a coding sequence\n",
    "        elif len(end_ind[0])==0:\n",
    "            non_coding_mutations[i]=1\n",
    "            original_aa.append('none')\n",
    "            new_aa_h.append('none')\n",
    "            new_aa_m.append('none')\n",
    "\n",
    "        #this leaves the coding mutations\n",
    "        else:\n",
    "            coding_mutations[i]=1\n",
    "\n",
    "            codon_idx = ind[0][0] #which codon is it\n",
    "            within_codon_idx = ind[1][0]#where does it fall in the codon (0,1,2)\n",
    "\n",
    "            #checking homology\n",
    "            #see if the codon falls in a region of alignment (checking list of codons/AAs that fall in region of hom.)\n",
    "            if codon_idx not in human_prot_align[gene_idx]:\n",
    "                non_homologous[i]=1\n",
    "                original_aa.append('none')\n",
    "                new_aa_h.append('none')\n",
    "                new_aa_m.append('none')\n",
    "\n",
    "            else: \n",
    "\n",
    "                aln_idx = human_prot_align[gene_idx].index(codon_idx)\n",
    "\n",
    "                human_aa_idx = human_prot_align[gene_idx][aln_idx] #finding index in mouse & human prot. sequence\n",
    "                mouse_aa_idx = mouse_prot_align[gene_idx][aln_idx]\n",
    "\n",
    "                #checking if flanking region is homologous\n",
    "                min_mouse = max(0, mouse_aa_idx-flank_size)\n",
    "                min_human = max(0, human_aa_idx-flank_size) #preventing weird errors\n",
    "                mouse_aa_flank = mouse_aa[gene_idx][min_mouse:mouse_aa_idx+flank_size+1]\n",
    "                human_aa_flank = human_aa[gene_idx][min_human: human_aa_idx+flank_size+1]\n",
    "\n",
    "                #flank matches = homologous at AA level\n",
    "                if mouse_aa_flank==human_aa_flank:\n",
    "                    #original_aa.append(human_aa[gene_idx][human_aa_idx]) #recording original amino acid\n",
    "                    homologous[i]=1\n",
    "\n",
    "                    ref_allele = mutation['Reference_Allele'].values[0]\n",
    "                    mut_allele = mutation['Tumor_Seq_Allele2'].values[0]\n",
    "\n",
    "                    h_codon = human_codon_seqs[gene_idx][human_aa_idx]\n",
    "                    m_codon = mouse_codon_seqs[gene_idx][mouse_aa_idx]\n",
    "\n",
    "                    #need to correct for strand differences with codon\n",
    "                    #checking human and mouse strand\n",
    "\n",
    "                    strand_h = df1[df1['gene']==gene_name]['strand'].values[0]\n",
    "                    if strand_h=='-': #if minus strand\n",
    "                        h_codon_true = h_codon.complement()\n",
    "                    elif strand_h=='+':\n",
    "                        h_codon_true = h_codon\n",
    "\n",
    "\n",
    "                    tx = homology_df[homology_df['gene']==gene_name]['mouse transcript'].values[0]\n",
    "                    strand_m = db_mouse[tx].strand\n",
    "                    if strand_m=='-': #if minus strand take complement\n",
    "                        m_codon_true = m_codon.complement()\n",
    "                    elif strand_m=='+':\n",
    "                        m_codon_true = m_codon\n",
    "\n",
    "\n",
    "                    aa_concordant, dna_concordant, ref_aa, mut_aa_h, mut_aa_m = mutation_modeling(h_codon, m_codon, codon_locs, codon_idx,human_aa_idx, mouse_aa_idx, within_codon_idx, mutation, gene_idx)\n",
    "\n",
    "\n",
    "                    aa_concordant_record[i] = aa_concordant\n",
    "                    dna_concordant_record[i] = dna_concordant\n",
    "                    original_aa.append(str(ref_aa))\n",
    "                    new_aa_h.append(str(mut_aa_h))\n",
    "                    new_aa_m.append(str(mut_aa_m))\n",
    "\n",
    "                    #if amino acid AND codon matches = homologous\n",
    "                    #if h_codon_true==m_codon_true:\n",
    "                    #    homologous[i]=1\n",
    "\n",
    "                    #if amino acid matches AND codon DOESN'T MATCH\n",
    "                    #else:\n",
    "                    #    aa_homologous_dna_non_homologous[i]=1\n",
    "\n",
    "\n",
    "                    #special case for deletion??\n",
    "                    #if mutation['Variant_Type'].values[0]=='DEL':\n",
    "                        #need to consider the entire size of the deletion (start and end site)????\n",
    "                    #NOT CURRENTLY CONSIDERING CASE WHERE INSERTION CAUSES DIFFERENT AA SEQUENCE???   \n",
    "\n",
    "                #flank doesn't match == non-homologous\n",
    "                else:\n",
    "                    non_homologous[i]=1\n",
    "                    original_aa.append('none')\n",
    "                    new_aa_h.append('none')\n",
    "                    new_aa_m.append('none')\n",
    "\n",
    "\n",
    "                    \n",
    "    df = pd.DataFrame(homologous, columns=['homologous'])\n",
    "    df['non_homologous']=non_homologous\n",
    "    df['coding_mutations']=coding_mutations\n",
    "    df['non_coding_mutations']=non_coding_mutations\n",
    "    df['excluded_mutations']=excluded_mutations\n",
    "    df['aa_concordant'] = aa_concordant_record\n",
    "    df['dna_concordant'] = dna_concordant_record\n",
    "\n",
    "    df['original_aa'] =  original_aa\n",
    "    df['new_aa_h'] = new_aa_h\n",
    "    df['new_aa_m'] = new_aa_m\n",
    "                  \n",
    "\n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Look at each mutation in dataset and quantifying homology\n",
    "- Iterating over different homology thresholds for region flanking the mutation of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/samgould/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:33: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "/Users/samgould/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:34: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "/Users/samgould/anaconda3/lib/python3.7/site-packages/Bio/Seq.py:2715: BiopythonWarning: Partial codon, len(sequence) not a multiple of three. Explicitly trim the sequence or add trailing N before translation. This may become an error in future.\n",
      "  BiopythonWarning)\n"
     ]
    }
   ],
   "source": [
    "for i in range(21):\n",
    "    flank_size=i\n",
    "    df = ortholog_PE(flank_size)\n",
    "    df['variant_type']=np.asarray(impact_data['Variant_Type'])\n",
    "    \n",
    "    path = '/Volumes/Sam_G_SSD/'\n",
    "    df.to_csv(path+'flanksize_' + str(i)+  '.csv')\n",
    "    \n",
    "    #these files are provided as well in the dropbox link"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "261.26837158203125px",
    "width": "551.635986328125px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
